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ABSTRACT 



Context. The shape of the stellar initial mass function (IMF) of a star cluster near its upper mass limit is a focal topic of investigation 
as it determines the high mass stellar content and hence the dynamics of the cluster at its embedded phase as well as during its young 
gas-free phase. The massive stellar content of a young cluster, however, can be substantially modified due to the dynamical ejections 
of the massive stars so that the present-day high-mass stellar mass function (hereafter MF) can be different than that with which the 
cluster is bora. 

Aims. In the present study, we provide a preliminary estimate of this evolution of the high-mass IMF of a young cluster due to early 
ejections of massive stars, using the Large Magellanic Cloud massive, young cluster R136 as an example. 

Methods. To that end, we utilize the results of the state-of-the-art calculations by Banerjee, Kroupa & Oh (2012) comprising direct 
N-body computations of realistic, binary-rich, mass-segregated models of R136. In particular, these calculations provide the ejection 
fraction of stars as a function of stellar mass. 

Results. We find that if the measured IMF of R136 is granted to be canonical, as observations indicate, then the "true" high-mass 
IMF of R136 at its birth must be at least moderately top-heavy when corrected for the dynamical escape of massive stars. 
Conclusions. The top-heaviness of the true high-mass IMF over the observationally determined one is a general feature of massive, 
young clusters where the dynamical ejection of massive stars is efficient. We discuss its implications and possible improvements over 
our current estimate. 

Key words. Galaxies: star clusters: general - Methods: numerical - Open clusters and associations: individual (R136) - Stars: 
kinematics and dynamics - Stars: luminosity function, mass function - Stars: massive 
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1. Introduction 



It is widely believed that stars form following an initial mass 
function (IMF) that has a universal form (Bast ian et al.| [2010 ; 
|Kroupa et al. 2012| l. In course of time however, the mass func- 
tion (MF) of a stellar population evolves as the stellar mem- 
bers are removed from the system. For bound stellar clusters, 
the low mass stars are naturally vulnerable to removal since 



they evaporate out of the system as the cluster relaxes ( Spitzer 



1987l |\£sperini & Heggie] [l997[ |Baumgardt & Makino[|20031 l 
The depletion of the low-mass stars from the bound system can 
be greatly augmented due to the expulsion of the residual gas 
that did not form stars, in the early phase of the cluster evolu- 
tion ( |Baumgardt & Kroupa[ |2007[ |Marks & Kroupaj |2010| l. A 
marked positive correlation of the global low-mass MF index 
(a) with the cluster's concentration has been pointed out for the 
first time by |De Marchi et al.| ( |2007| i. |Baumgardt et al.| (2008 ) 
showed that such a correlation can be explained by clusters born 
with a canonical IMF which are primordially mass-segregated 
but are filling their tidal radii. [Marks & Kroupa ( 2010[ ), how- 
ever, demonstrated that the De Marchi trend can as well be ex- 
plained by residual gas expulsion from initially compact and 
mass-segregated clusters without introducing a variable low- 
mass IMF. This result is consistent with the compactness of 
freshly formed clusters (Ma rks & Kroupa] |2012| l, whereas the 



tidal-radius -filling ansatz would require the clusters to condense 
from a cloud with large radius matching the tidal radius. 

On the other hand, the distribution of the most massive stars 
bound to a cluster is influenced by the ejection of the massive 
members via dynamical encounters (Clarke & Pringle 1992] >. 
Massive stars in a young cluster, which are usually in bina- 
ries, segregate in the cluster's center, where their density can be 
greatly augmented. The massive stars therefore frequently in- 
teract through super-elastic binary-single and binary-binary en- 
counters to eject each other from the cluster, that way deplet- 
ing the massive stellar population as the cluster evolves. In gen- 
eral, more massive stars are more centrally concentrated and 



hence are more severely depleted (Pflamm-Altenburg & Kroupa 
2006). This is found to be true in theoretical studies of mas 



sive r unaway sta rs from star clusters by seve ral authors (Fujii & 
Portegies Zwart 2011 Banerjee et al.||2012a) . 



The increasing depletion of (bound) massive stars with in- 
creasing stellar mass in star clusters has remarkable implications 
on its IMF. If a young cluster is found with the massive end of 
its MF in the canonical form (Kroupa |2001| l at its present age, 
then the upper end of its IMF has to be top-heavy to achieve 
the canonical form at the current age (taking into account the 
mass loss due to stellar evolution; see Sec.[2]i. In general, given 
a measured MF of a young cluster near the massive limit, the 
corresponding stellar mass distribution at birth, i.e., the actual 
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Fig. 1. The fraction f e] (M) of runaway stars for M > 10M G at t x 3 Myr cluster age as obtained from the four computations of 
Paper I (c.f. Fig. 4 of Paper I). We use the same symbols as in Fig. 4 of Paper I to indicate that they are obtained from direct N-body 
computations starting with different initial realizations of the R136 model employed in Paper I. The central solid line represents the 
best-fit line to the data that asymptotes to unity (see Sec. [4}. The upper and the lower dashed lines limit the error in the fitting. The 
dotted line is the best-fit single power-law to the entire data which is truncated at f e j(M) - 1 (see Sec. |4jl. 



IMF, has to be "top-heavier" depending on the age of the cluster 
due to stellar-dynamical ejections. 

In the present work we provide a preliminary quantification 
of the top-heaviness of the IMF of a star cluster as inferred from 
the above considerations. To that end we utilize the ejected frac- 
tion of massive stars, / e j(M), as a function of the mass M of 
the ejected stars from model clusters as computed by Banerjee 
et al. (2012a; hereafter Paper I). The above quantity is defined 
as the ratio of the number of stars, within a bin around mass 
M, moving away from the cluster beyond R > 10 pc from the 
cluster's center to the total number of stars in the whole system, 
i.e., including all the bound and the ejected members, within the 
same mass bin. Note that f ej {M) was denoted by g(M) in Paper 
I but we stick to this new symbol in this paper for better clarity. 
The above authors studied the runaway O-stars from the Large 
Megallanic Cloud (LMC) cluster R136 by computing the dy- 
namical evolution of realistically modelled R136-type clusters 
using direct N-body integration. We point out the essentials of 
these computed models of Paper I in Sec. [3] 



pointed out in Sec. [T £(M) needs to be further corrected for the 
stars removed from the bound cluster to determine the true IMF 
of the cluster, ^ me (M), which is the bound stellar mass function 
of the cluster at its birth. To correct ^(M) over the entire mass 
range, one should, in principle, reassemble the cluster including 
all the stars that escaped from the cluster. If the cluster is formed 
mass-segregated (as considered in Paper I and inferred to be true 
for several young clusters; see, e.g., |Bonnell & Dav ies 1998), 



then the massive stars reside in the cluster's center from the very 
beginning of its secular evolution and they can only be ejected 
by dynamical encounters {i.e., they cannot be removed by the 
residual gas expulsion or the external galactic field). Therefore, 
^ime(M) can be obtained from %{M) by only correcting for the 
runaway stars, for the massive end. 

In Paper I, the ejection fraction / e j(M, t) is found to increase 
significantly with stellar mass M for M > 1OM at the cluster 
age of t « 3 Myr (c.f. Fig. 4 of Paper I). The number of ejected 
stars over the mass range M to M + dM is, 



n ej (M, t)dM = / ej (M, f)^ me (M)dM. 



(1) 



2. Dynamical evolution of the bound stellar mass 
function of a star cluster 

In reality, the present-day observationally determined bound 
stellar mass function of a young star cluster is corrected for un- 



resolved multiple companions (Kroupa, Gilmore & Tout 1991 



|Weidner, Kroupa, & Mas chberger, 2009J, and for the wind mass 
loss of the stars to determine the so-called cluster IMF, £(M). 
The correction for unresolved multiple components is negligible 
for massive stars (Weidner, Kroupa, & Maschberger 2009}. As 



Here, the decrement of the zero-age mass, AM(M, t), is not in- 
cluded since in evaluating / e j(M, t), the stellar masses at age t 
are used for both the total population between M to M + dM and 
the corresponding ejected population so that the same mass shift 
(AM) is applicable to all the three functions n e j, / e j and £ true . For 
brevity, we also suppress the explicit t dependence in Eqn. ([T]). 

The number of stars in the mass range M to M + dM remain- 
ing bound to the cluster is therefore, 



{ lme (M)dM - n e] (M)dM = £ rae (M){l - f e] (M)}dM, 



(2) 
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Fig. 2. The "true" IMF, ^ a -u e {M), as corrected for the ejected stars assuming the "observed" IMF, £(M), to be the Salpeter/canonical 
power law %(M) oc M~ 23 (represented by the dash-dotted line). The true IMF is related to the observed IMF by Eqn.[4] The upper 
dashed, the lower dashed, the central solid and the dotted lines correspond to those in Fig.[T]respectively. It can be noted that the 
correction of the observed IMF due to the dynamical ejection of massive stars can amount to a moderately top-heavy true IMF (the 
dashed and the solid lines). However, the data in Fig.[T]is also supportive of a substantially higher top-heaviness (the dotted line). 
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Fig. 3. The ratio of the "true" IMFs, ftmeW. in Fig. [2] to the "observed" IMF, £(M), the latter being assumed to be a 
Salpeter/canonical power law £(M) oc M~ 21 . The upper and the lower dashed, the central solid and the dotted lines correspond 
to those in Fig.[T]respectively. 
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where we have used Eqn. ([TJ. The right-hand-side of Eqn. Q is 
the "observed" (high-mass) IMF of the cluster without including 
the ejected stars (but corrected for stellar evolution). In other 
words, 



£(M)=£ rae (M){l-/ ej (M)}, 
or, 

£(M) 



ftme(M) 



I -MM)' 



(3) 



(4) 



Eqn. Q connects the observed cluster IMF, ^(M), with its true 
IMF, ftueW 



3. Previous computations 

State-of-the-art calculations have been performed in Paper I to 
evolve model star clusters that mimic R136 in terms of observ- 
able parameters. The primary objective of these calculations, 
which are so far the largest sized ones with fully mass-segregated 
hard massive binaries (see below), was to study the ejected OB- 
stars from R136, in particular, whether massive runaway stars 
like VFTS 682 and 30 Dor 016 can indeed be ejected dynami- 



cally from R136 as suspected (|Bestenlehner et al. 2011 Evans 



|et al.||2010] ). These computations demonstrated good agreement 
of the kinematics of the massive ejecta with those of these noted 
runaways. Furthermore, they suggest that the tight massive bi- 
naries that drive these runaways can merge due to the dynam- 
ical interactions leading to the formation of "super-canonical" 



(Kroupa et al. 2012| l stars, i.e., stars with masses exceeding the 
canonical upper limit of 15OM (see Paper I and Bane rjee et] 
al.|2012b I, in agreement with the observed stellar population in 
R136 ( jCro wther et aLj |2010| ). In this study, we continue to utilize 
the results of these computations. The initial conditions and the 
method of these calculations are accounted in detail in Sec. 2 of 
Paper I which we briefly recapitulate here. 

The above computations comprise direct N-body integra- 
tions of Plummer spheres ( Heggie & Hut 2003 |Kroupa| [2008 ) 
with parameters conforming with those of R136. The initial mass 
of the Plummer spheres is M e ;(0) « 1O 5 M which is an upper 
limit for R136 (Crowth er et al.| [2010 > and the half mass radii 
are taken to be r/,(0) » 0.8 pc. The clusters are made of stars 
withZAMS masses drawn from a canonical IMF (Kroupa, 2001 ) 
over the range 0.08M G < m s < 15OM and are of metallicity 
Z = 0.5Z G . With the mean stellar mass of m Q.6M Q , each cluster 
comprises w 1.7 x 10 5 stars. 

As for the primordial binary population, stars with m s > 
5M Q are all in binaries while all lighter stars are kept single. 
This termination of the binary population is due to computa- 
tional ease; binaries bottleneck the calculation speed of direct 
N-body integration significantly so that adapting a full spectrum 
of primordial binaries (i.e., 100% initial binary fraction; Marks 
& Kroupa|2012| l in models of the size that we compute becomes 
prohibitive. Disregarding the binary population for m s < 5M of 
course does not affect the ejection of massive stars and merg- 
ers of massive binaries significantly (see Paper I for details). 
Following the observed period distribution of O-star binaries 
(Sana & Evans 2011 1, the orbital periods of the binaries having 
primary masses m s > 2OM are chosen from a uniform distribu- 
tion between 0.5 < log 10 P < 4, where P is the orbital period in 
days. For primary masses 5M < m s < 2QM Q , the "birth period 



These distributions, however, are not modified via "eigenevolu- 
tion" (Kroupa 1995) as the interactions between massive pre- 
main-sequence stars are not yet quantified. The initial configura- 
tions are fully mass-segregated (Baumgardt et alT| |2QQ8) > so that 
the massive binaries concentrate within the clusters' central re- 
gion mimicking primordial mass segregation which is inferred 
for several Galactic globular clusters and young star clusters. 

To perform the direct N-body computations, the state-of-the- 
art "NBODY6" integrator ( Aarseth"l |2003| l has been used. In ad- 
dition to integrating the particle orbits using the highly accu- 
rate fourth-order Hermite scheme and dealing with the diverg- 
ing gravitational forces in close encounters through regulariza- 
tions, NBODY6 also employs the analytical stellar and binary 
evolution recipes of Hurley et al. (2000 2002). The models are 
computed for m 3 Myr four times, each with a different random 
number seed initialization. 



4. Ejected fraction of massive stars / e j(M): a 
top-heavy IMF in R136 

Figure 4 of Paper I shows the ejected fraction of stars as a func- 
tion of stellar mass M, where a significant and monotonic in- 
crease of fej(M) can be seen beyond M > 1OM . As explained 
in Paper I, this trend is a direct consequence of mass segrega- 
tion; more massive systems are more concentrated towards the 
cluster's center and hence interact more efficiently (see below). 
Fig.[T|re-plots the f e j(M) at the cluster age t w 3 Myr, as obtained 
from the computations of Paper I, for M > 1OM where we use 
the same symbols as in Figure 4 of Paper I to differentiate among 
the data from different computed models. 

A possible form of the fitting function to the data in Fig. [T] 
would be asymptotic to 1 with increasing M as / e j(M) < 1 by 
definition but monotonically increasing as the trend of the data 
suggests. Due to the scatter of the points in Fig. [T] over a large 
range and to ensure reliability, we do a least-squares fit to the 
logarithm of the data in Fig. [T] We adopt the exponential func- 
tional form, 



log(/ ej (M)) = -Be P x (B >0,p>0), 



(5) 



to best-fit with the logarithm of the data-coordinates of Fig. [T] 
This functional form asymptotically increases to 0, i.e., ap- 
proaches 1 when transformed into the linear scale. We use all 
the 28 data-points of the above figure to obtain the best fit val- 
ues B = 8.03 + 1.42 and p = 1.20 + 0.13. The central solid 
line in Fig. [TJrepresents the best-fit line to the data and the two 
dashed lines represent the upper and the lower bounds of the fit 
corresponding to the errors (1-cr or the standard errors) in B and 
P- 

As such, the data-points in Fig.[T]can be found to follow an 
overall power-law index of m 1.5 if one fits a power-law to them. 
A basic physical argument that supports such an increment is as 
follows. The massive single stars get ejected by super-elastically 
(Heg gie & Hut| |2003[ ) encountering with binaries within the 
cluster (as noted in Paper I, most of the ejected entities are sin- 
gle stars, given the initial distribution of the binaries' orbital pe- 
riods adapted her^J. These single stars are liberated from bina- 
ries either by binary-binary encounters or by direct disruption of 
binaries and they can, in general, be expected to remain mass- 
segregated within the cluster. A given segregated group of stars 



distribution" of Kroupa ( 1995 1 is adapted. The eccentricity dis- 
tribution is chosen to be thermal ( |Spitzer| [1987| [Krou pa 2008 ). 



1 Note that ejected single stars may in reality be binaries if the mod- 
els would incorporate a large fraction of O-stars being in very tight (< 
few au) binaries within wider triple or quadruple systems with outer 
orbital periods as chosen here (see Chini et al. 2012). 
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of mass M remain confined within a radius r oc M~ 1 ' 2 of the clus- 
ter (see Chapter 16 of Heggi e~& Hut|2003"] ) so that their density 
oc r~ 3 oc M 3/2 . This density is also proportional to the single- 
stars' encounter rate with the binaries and hence to the proba- 
bility of the formers' ejection from the cluster. Hence, one can 
expect that the ejection fraction f e j(M) oc M 1 ^ 2 which is indeed 
close to the best-fitted power-law index as obtained here. 

The dotted line in Fig. [T] represents the best-fit power law 
/ e j(M) = AM k (truncated at / e j(M) = 1). By doing a linear least- 
squares fit to the logarithm of the data-points, we obtain the best- 
fit values A = i(r 3 - 65±() - 22 and A; = 1.46+0.13. Unfortunately, due 
to the sparsity of data-points near the massive end it is presently 
difficult to determine the true functional behavior of f e j(M) over 
the entire data of Fig.[T] A larger number of such computations 
would help to settle this issue. 

The true IMFs, £ tlue (M), corresponding to the f e j(M) lines 
in Fig. [T] as obtained from Eqn. (j4j), is shown in Fig. [2] where 
the observed IMF is taken to be the canonical one, i.e., follow- 
ing the Salpeter form £(M) oc M~ 23 for M > 10M G (shown 
by the dash-dotted line). The upper dashed, central solid, lower 
dashed and the dotted lines in Fig.|2]correspond to those in Fig.[T] 
respectively. Noticeably, the IMF of R136 is indeed observed 
to be very close to canonical over both low-mass (Andersen et 
|al.|2009| 1.1 - 10M G ) and high-mass ( |Massey & Hunter|1998| 
2.8 - 120M G ) ranges. 

It can be seen in Fig. [2] that the true IMF, ^tme(M), is top- 
heavy near its massive end even for M < 150M Q . In other words, 
the computed ejection fraction of massive stars from R136-type 
model clusters (in Paper I) indicates from a moderate to a sub- 
stantial top-heaviness in the high-mass IMF of R136. This point 
becomes more vivid in Fig. [3] where we plot the ratio of the true 
IMF, ^true(M), to the Salpeter law £(M) oc AT 23 . 

Notably, the fitted bounds to the data-points in Fig. [T] (the 
dashed lines) represent the effect of statistical variation of the 
ejection-fraction, / e j(M), between the computed cluster models. 
Such cluster to cluster variation affects the top-heaviness of the 
inferred true IMF as can be seen in the corresponding bounds of 
^tmeiM) m Fig. [2] The inherent statistical fluctuations can result 
in the inference of a marginally top-heavy to a moderately top- 
heavy true IMF (c.f. Fig. |3j. In other words, if the present-day 
observed IMF in young clusters is canonical, then there has to 
be a cluster to cluster stochastic variation of the top-heaviness 
of the true IMF with which they are born. Further detailed study 
(c.f. Sec. [5]) is needed to quantify this better. [Dib et al. (2010) 
has studied the growth of massive stellar cores by gas accretion 
and the subsequent formation of stars including the effects of 
stellar wind feedback on the protocluster which indicates cluster- 
to-cluster variation of the stellar IMF (see also |Dib et al.|2007) . 



5. Discussions and outlook 

The above inference of top-heaviness near the upper end of the 
"true" IMF is not limited to R136 but is true for any massive 
cluster that is observed to have a normal (i.e., canonical) IMF. 
The extent of the top-heaviness of the true IMF for a particular 
cluster, of course, depends on its ejection fraction / e j(M). It is to 
be noted that the above analysis is only of the lowest order and is 
essentially a qualitative inference. This is because had we started 
the N-body computations with the inferred top-heavy true IMF 
as corrected for the runaways, f e j(M) would have been larger 
than the presently used one which is obtained from a canonical 
IMF. Therefore, in general, the (top-heavier) true massive end 
of the IMF can be determined quantitatively from an observed 



IMF by iteratively applying the correction in Sec. |4] for the dy- 
namical ejection of the massive stars which would tend to aug- 
ment the top-heaviness. Such a study would reveal relationships 
between the top-heaviness of a cluster's IMF with its initial pa- 
rameters such as the total mass, the degree of primordial mass 
segregation and the concentration subject to the constraint that a 
nearly canonical mass distribution is reached at the present day 
which is true for many Galactic young clusters. Such relation- 
ships would, in turn, have fundamental implications of massive 
star formation in a clustered environment (Bon nell et al.| [l998; 
|Zinnecker & Yorke||2007) . 

It is important to note that the top-heaviness that has been 
inferred above for the massive (true) IMF is influenced by 
the adopted initial complete mass-segregation (see Sec. [3) that 
provides the maximum efficiency in ejecting massive stars at 
young ages. However, such an assumption of primordial mass- 
segregation is in accordance with the observations of young star 
clusters; see Paper I and references therein (and, e.g., Bonnell 
|& Davies]|1998| l. As pointed out above, an interesting outlook 
would be to study the degree of top-heaviness of the true IMF 
as a function of the initial degree of segregation. The timescale 
of segregation of the massive stars/binaries is also shortened by 
the cluster's concentration. Therefore, it would also be worth- 
while to study the dependence of the effect on the cluster's ini- 
tial concentration. A physical basis of the true IMF being highly 
top-heavy can be found in |Dib et aL] ( |2010| l. In this study of 
the appearance of the stellar IMF via accretion onto proto-stellar 
cores including wind feedback, these authors do find a substan- 
tial flattening over the massive end of the IMF. 

In Fig. |TJ the best-fit / e j(M) plot extends beyond the canoni- 
cal upper limit of 150M o . As explained in Paper I and Banerjee 



|et al.| ( |2012b| l, this is to take into account the "super-canonical' 
stars with zero-age masses between 150 - 300M G that form 
via mergers of massive binaries due to dynamical encounters. 
Although the corrected IMF, £ tme (M), is plotted in Fig.|2]for M > 
15OM consistently with Fig.|T] only % tme (M) for M < 15OM is 
relevant as the cluster's IMF. 

In conclusion, our above analysis with the example of R 136 
indicates that taking into account the dynamical ejection of mas- 
sive stars, a young cluster with a seemingly normal or canonical 
IMF can actually have formed with an IMF which is at least 
moderately top-heavy near its high-mass end. This may well 
constitute the first ever detection of the flattening of the IMF 
as predicted by |Dib et al.| ( |2010"| ). For completeness, we mention 
that by tracing the present-day constitution of globular clusters 
to their initial configurations, a trend with birth density and the 
power-law index of the IMF has emerged such that the IMF is in- 
creasingly top-heavy with increasing star-forming dens ity above 
about > 10 5 M o pc~ 3 on a cluster-forming spatial scale ( |Marks et] 



al. 2012[). The present results thus add to the increasingly solid 



evidence that the IMF becomes top-heavy in star bursts. 
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